
      subroutine makepic(l,n,xpts,ypts1,ypts2,ypts3,xb,xe,yb,ye)
      character*11 string(7)
      dimension xpts(n),ypts1(n),ypts2(n),ypts3(n),x(200),y(200),
     >   ybar(200)
      data string/'10!SP!1!RE','10!SP!0!RE','10!SP!-1!RE',
     >   '10!SP!-2!RE','10!SP!-3!RE','10!SP!-4!RE','10!SP!-5!RE'/
      
      call pspace(xb,xe,yb,ye)
      call hrdchr(0)
      call mskchr(1)
      i = 1
      ymax = 0.0
      ymin = 1e10
      do j = i, n
         if (ypts3(j).gt.ymax) ymax = ypts3(j)
         if (ypts2(j).gt.ymax) ymax = ypts2(j)
         if (ypts1(j).gt.ymax) ymax = ypts1(j)
         if (ypts3(j).lt.ymin) ymin = ypts3(j)
         if (ypts2(j).lt.ymin) ymin = ypts2(j)
         if (ypts1(j).lt.ymin) ymin = ypts1(j)
      end do 
      ymax = ymax * 1.1
      ymin = ymin * 0.8
      xmin = 0.0
      xmax = 180.0
      if (l.le.2) then
         call mapyl(xmin,xmax,ymin,ymax)
         call mask(xmin-40.0,xmin-1.0,ymin/1.2,ymax*2.0)
         call scayli(20.0)
         call unmask(0)
         do j = 1, 7
            r = 10.0**(2-j)
            if (r.lt.ymax.and.r.gt.ymin)
     >         call plotst(xmin-12.0,r,string(j))
         end do 
      else
         if (l.eq.4) ymin = -0.2
         call map(xmin,xmax,ymin,ymax)
         call scalsi(20.0,0.2)
      end if 
      call curveo(xpts,ypts3,i,n)
      call broken(10,5,10,5)
      call curveo(xpts,ypts2,i,n)
      call broken(5,5,5,5)
      call curveo(xpts,ypts1,i,n)
C  set lines to be full again
      call full
      if (l.eq.1) then
         call rexp('exp.2s',x,y,ybar,npts,1.0)
         call erry(x,y,ybar,1,npts,229)
      else if (l.eq.2) then
         call rexp('exp.2p',x,y,ybar,npts,1.0)
         call erry(x,y,ybar,1,npts,229)
      else if (l.eq.3) then         
         call rexp('lr.exp.jw',x,y,ybar,npts,1.0)
         call erry(x,y,ybar,1,npts,229)
         call rexp('lr.exp.ew',x,y,ybar,npts,1.0)
         call erry(x,y,ybar,1,npts,227)
      else 
         call rexp2('lr.exp.jw',x,y,ybar,npts,1.0)
         call erry(x,y,ybar,1,npts,229)
         call rexp2('lr.exp.ew',x,y,ybar,npts,1.0)
         call erry(x,y,ybar,1,npts,227)
      end if 
      return
      end
      
      subroutine rexp(file,x,y,ybar,n,c)
      character file*(*)
      dimension x(200),y(200),ybar(200)
      
      open(99,file=file)
      i = 1
 15   read(99,*,end=25) x(i),y(i),ybar(i)
      y(i) = y(i) * c
      ybar(i) = ybar(i) * c
      i=i+1
      go to 15
 25   n=i-1
      close(99)
      return
      end
      
      subroutine rexp2(file,x,y,ybar,n,c)
      character file*(*)
      dimension x(200),y(200),ybar(200)
      
      open(99,file=file)
      i = 1
 15   read(99,*,end=25) x(i),a,b,y(i),ybar(i)
      y(i) = y(i) * c
      ybar(i) = ybar(i) * c
      i=i+1
      go to 15
 25   n=i-1
      close(99)
      return
      end
      
      subroutine rdata(file,x,y,n)
      character file*(*)
      dimension x(200),y(200)
      open(99,file=file)
      i = 1
 15   read(99,*,end=25) x(i),y(i)
      i=i+1
      go to 15
 25   n=i-1
      close(99)
      return
      end
            
C  The following routines generate error bars with open symbols.
C  Written by Igor Bray (1992)
      Subroutine erry(x,y,ybar,istart,istop,icode)

C  uses GHOST routines to plot vertical error bars on a graph.
C  If icode is zero or greater, a marker symbol is plotted at
C  each x,y point.

      Dimension x(istop),y(istop),ybar(istop)
      Dimension xtemp(2),ytemp(2)

      if(icode.ge.0) call ptplot(x,y,istart,istop,icode)
      do 10 i=istart,istop
         xtemp(1)=x(i)
         xtemp(2)=x(i)
         ytemp(1)=y(i)+ybar(i)
         ytemp(2)=y(i)-ybar(i)
         call getminmax(x(i),y(i),xmin,xmax,xminp,ymin,ymax,yminp)
         call mask(xmin,xmax,ymin,ymax)
         if (ybar(i).gt.y(i)-yminp) call ptgraf(xtemp,ytemp,1,2,45)
         call unmask(0)
10    continue
      return
      end

      Subroutine errx(x,y,xbar,istart,istop,icode)

C  uses GHOST routines to plot horizontal error bars on a graph.
C  If icode is zero or greater, a marker symbol is plotted at
C  each x,y point.

      Dimension x(istop),y(istop),xbar(istop)
      Dimension xtemp(2),ytemp(2)

      if(icode.ge.0) call ptplot(x,y,istart,istop,icode)
      do 10 i=istart,istop
         ytemp(1)=y(i)
         ytemp(2)=y(i)
         xtemp(1)=x(i)+xbar(i)
         xtemp(2)=x(i)-xbar(i)
         call getminmax(x(i),y(i),xmin,xmax,xminp,ymin,ymax,yminp)
         call mask(xmin,xmax,ymin,ymax)
         if (xbar(i).gt.x(i)-xminp) call ptgraf(xtemp,ytemp,1,2,192)
         call unmask(0)
10    continue
      return
      end

      Subroutine errxy(x,y,xbar,ybar,istart,istop,icode)
C  uses GHOST routines to plot both vertical and horizontal error
C  bars on a graph. If icode is zero or greater, a marker symbol
C  is plotted at each x,y point.

      Dimension x(istop),y(istop),ybar(istop),xbar(istop)
      Dimension xtemp(2),ytemp(2)

      if(icode.ge.0) call ptplot(x,y,istart,istop,icode)
      do 10 i=istart,istop
         xtemp(1)=x(i)
         xtemp(2)=x(i)
         ytemp(1)=y(i)+ybar(i)
         ytemp(2)=y(i)-ybar(i)
         call getminmax(x(i),y(i),xmin,xmax,xminp,ymin,ymax,yminp)
         call mask(xmin,xmax,ymin,ymax)
         if (ybar(i).gt.y(i)-yminp) call ptgraf(xtemp,ytemp,1,2,45)
         ytemp(1)=y(i)
         ytemp(2)=y(i)
         xtemp(1)=x(i)+xbar(i)
         xtemp(2)=x(i)-xbar(i)
         if (xbar(i).gt.x(i)-xminp) call ptgraf(xtemp,ytemp,1,2,192)
         call unmask(0)
10    continue
      return
      end

      subroutine getminmax(x,y,xmin,xmax,xminp,ymin,ymax,yminp)
      dimension idata(10),rdata(12)

C  get the magnification of the character
      call enqchr(idata,rdata)
      mag = idata(2)
C  The size in ND space of half the character
      size = mag * 0.001 / 7.0 * 4.0 / 2.0
C  Decide on the min span of error bars in order for them to be printed
      sizep = size * 1.5
C  Get the mappings
      call enqmap(idata,rdata)

      if (idata(1).le.2) then
         const = (rdata(2) - rdata(1)) / (rdata(10) - rdata(9)) 
         xmin = x - const * size 
         xminp = x - const * sizep
         xmax = x + const * size 
      else
         const = (log(rdata(2)) - log(rdata(1)))/(rdata(10) - rdata(9))
         xmin = x / exp(const * size)
         xminp = x / exp(const * sizep)
         xmax = x * exp(const * size)
      end if

      if (idata(1).eq.1.or.idata(1).eq.3) then
         const = (rdata(4) - rdata(3)) / (rdata(12) - rdata(11)) 
         ymin = y - const * size 
         yminp = y - const * sizep
         ymax = y + const * size 
      else
         const = (log(rdata(4)) - log(rdata(3)))/(rdata(12) - rdata(11))
         ymin = y / exp(const * size)
         yminp = y / exp(const * sizep)
         ymax = y * exp(const * size)
      end if
      end
      
